clear all;
%% Code explanation
% 170305 DZ 
% 190709 DZ updated to calculate 50%Rn and 1%Rn
% 190712 DZ offset for resistance

%% Define the resistance to temperature conversion equation
T_RuOx=@(R)1./exp(-1.001200439904+1.175279377873*log(R-2.2)+0.021010648114*(log(R-2.2)).^2-0.024007137972*(log(R-2.2)).^3+0.002054398372*(log(R-2.2)).^4);
%
%% Set the parameter

Boffset=0.06;
L_W=24.2/53.3;
ratio=0.5;

%% Read the data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% T=75;
%  T=[10 25 37 50 62 75 87 100 112 125 137 150 175 200]';

root=('D:\Documents2\Rawdata\Sn\171219 MPI mK\20180330 2Sn15PbTe\'); %

filenames={'20180223_0uW';'20180224_10uW';'20180225_20uW';'20180226_40uW';...
    '20180226_70uW';'20180228_110uW';'20180228_140uW';'20180302_170uW';...
    '20180304_210uW';'20180306_250uW'; '20180327_270uW';'20180307_300uW';'20180326_380uW';'20180322_400uW';'20180326_470uW';'20180323_550uW';...
    '20180324_700uW';'20180324_900uW';'20180325_1100uW';'20180325_1300uW'};


 Rnorm=6.4284; % normal state resistances of the two samples

NT=length(filenames);


colorset=jet(NT);
figure

Tmc=zeros(NT,1);
Tmc_err=zeros(NT,1);
Hc2=zeros(NT,1);
Rc=zeros(NT,1);
% Tmax1=0.92;
% Tmax2=1.05;
% Tmc(NT+1,1)=Tmax1;
% Tmc(NT+1,2)=Tmax2;
%%
kk=NT;
 file1=[root,filenames{kk},'_wroots_inplane_waitmeas.dat'];
        fid=fopen(file1);
        data1 = textscan(fid,'%s %s %s %s','commentStyle','#');
        data1=data1';
        for ii=1:size(data1,1)
            data1{ii,:}=str2double(data1{ii,:});
        end       
        fclose(fid);
    
     
          B=round(data1{1}*10)/10;
          R1=(data1{2}/L_W+0.05)/50;
         Bset=(min(B):0.1:max(B))';
         Bset=round(Bset*10)/10;
         RT1=zeros(length(Bset),1); % T and R1
   
         for ii=1:length(Bset) % for every fixed B field, find the minimal R     
            [RT1(ii),j1]=min(R1(B==Bset(ii)));
         end
         
          Bset=Bset-Boffset;      
          indexB=find(Bset>=1.4);
          p1=polyfit(Bset(indexB),RT1(indexB),1);
 
          subplot(2,2,2),plot(Bset,ratio*(p1(1)*Bset+p1(2)),'-r');
          hold on
%%
for kk=1:NT
   
    file1=[root,filenames{kk},'_wroots_inplane_waitmeas.dat'];
        fid=fopen(file1);
        data1 = textscan(fid,'%s %s %s %s','commentStyle','#');
        data1=data1';
        for ii=1:size(data1,1)
            data1{ii,:}=str2double(data1{ii,:});
        end       
        fclose(fid);
    
     
          B=round(data1{1}*10)/10;
          R1=(data1{2}/L_W+0.05)/50;
          
          T=T_RuOx(data1{3}*100);
%           theta=data1{4};
            Tmc_err(kk,1)=3*std(T);
            
         B0=(min(B):0.002:max(B))';
         Bset=(min(B):0.1:max(B))';
         Bset=round(Bset*10)/10;
         RT1=zeros(length(Bset),2); % T and R1
      
    
         
         for ii=1:length(Bset) % for every fixed B field, find the minimal R
            indexB=find(B==Bset(ii));
            [RT1(ii,2),j1]=min(R1(indexB));
            RT1(ii,1)=T(indexB(j1));     
         end
          RRT1=interp1(Bset,RT1(:,2),B0,'spline');
          TT1=interp1(Bset,RT1(:,1),B0,'linear');
          Bset=Bset-Boffset;
          B0=B0-Boffset;
           subplot(2,2,1),plot(Bset,RT1(:,2),'s','Color',colorset(kk,:)','MarkerSize',3);
          hold on
            subplot(2,2,2),plot(B0,RRT1,'-','Color',colorset(kk,:)','LineWidth',1);
          hold on
          


          index1=find((RRT1-ratio*p1(1)*B0-ratio*p1(2)<=0));
          
         if isempty(index1)==0
             Hc2(kk)=B0(index1(length(index1)));           
          Tmc(kk)=TT1(index1(length(index1)),1);
          Rc(kk)=RRT1(index1(length(index1)));
         else
             Hc2(kk)=NaN;
             Tmc(kk)=NaN;
             Rc(kk)=NaN;
         end
          title={'B(T)','R(kOhm)','T(K)','averaged T(K)'};
        AA=[Bset RT1(:,2) RT1(:,1) mean(RT1(:,1))*ones(length(Bset),1)];
       
         xlswrite('FigS6B.xlsx',AA,kk);
       xlswrite('FigS6B.xlsx',title,kk);
          
         
%     
end



figure
plot([Tmc' 0.3795]',[Hc2' 0]','or');
hold on


